Summary of experiments DTU

#code for DTU analysis can be found in docs DTU.Rmd. Only follow up analyses are showed in html.

empirical FDR < 0.05

rm(sumExp)
# summarise

# empirical FDR < 0.05
map_df(all_res, ~{group_by(.x, FDR = empirical_FDR < 0.05) %>% tally() %>% filter(FDR == TRUE)  }, .id = "contrast")

regular FDR < 0.05

# regular FDR < 0.05
map_df(all_res, ~{group_by(.x, FDR = regular_FDR < 0.05) %>% tally() %>% filter(FDR == TRUE) }, .id = "contrast")

Upsetplot DTU across conditions

dtu_genes <- map(all_res, ~{ .x %>% filter(empirical_FDR < 0.05) %>% pull(gene) })

#UpSetR::upset(fromList(dtu_genes),nintersects = NA)
UpSetR::upset(fromList(dtu_genes),nintersects = NA,nsets = 7)

sig_res <- map_df(all_res, ~{.x %>% filter(empirical_FDR < 0.05, abs(estimates) > 1)}, .id = "contrast")

age_res <- filter(all_res$age, empirical_FDR < 0.05)

stimulation_res <- map_df(all_res, ~{.x %>% filter(empirical_FDR < 0.05, abs(estimates) > 1)}, .id = "contrast" ) %>% filter(contrast != "age")

Volcano plots stimulation

all_res %>%
  map_df( ~{.x %>% filter(pval < 0.05)}, .id = "contrast") %>%
  filter(!contrast %in% "age") %>%
  ggplot(aes( x = estimates, y = -log10(pval), colour = empirical_FDR < 0.05 )) +
  geom_point(size = 1) + scale_colour_manual(values = c("black", "red")) + 
  facet_wrap(~contrast, scales = "free") +
  xlim(-4,4) +
  theme_bw()
## Warning: Removed 4 rows containing missing values (geom_point).

Volcano plot age

# age only
all_res %>%
  map_df( ~{.x %>% filter(pval < 0.05)}, .id = "contrast") %>%
  filter(contrast == "age") %>%
  ggplot(aes( x = estimates, y = -log10(pval), colour = empirical_FDR < 0.5 )) +
  geom_point(size = 1) + scale_colour_manual(values = c("black", "red")) + 
  facet_wrap(~contrast, scales = "free") +
  xlim(-0.1,0.1) +
  theme_bw()
## Warning: Removed 6727 rows containing missing values (geom_point).

Boxplot for TYROBP in LPS

markers = "~/Documents/MiGASti/docs/2nd_pass/TYROBP.xlsx"
markers = read_excel(markers, col_names = TRUE) 
markers = as.data.frame(markers)

isoform_id <- rownames(tx_ratio)
copy <- cbind (isoform_id, tx_ratio)
ratio <- merge(copy, markers, by = "isoform_id")
sample <- colnames(ratio)
rownames(ratio) = ratio$isoform_id
metadata <- read.table("~/Documents/MiGASti/Databases/metadata.txt")
#set rownames to Sample
row.names(metadata) <- metadata$Sample 
setwd("~/Documents/MiGASti/Databases")
exclude <- read.table("samples2remove.txt")
exclude <- exclude$x
ratio = ratio[, !colnames(ratio) %in% exclude] 
#Excludes the samples from filters. 
#dim(genes_tpm_filt)
metadata_filt = metadata[ !(rownames(metadata) %in% exclude), ]
length(metadata_filt)
## [1] 38
samples_LPS = metadata_filt$Sample[metadata_filt$Stimulation %in% "LPS"]  
samples_unstim = metadata_filt$Sample[metadata_filt$Stimulation %in% "unstim"] 

metadata4deg = metadata_filt[rownames(metadata_filt) %in% c(samples_LPS,samples_unstim) ,]

ratio_LPS_unstim = as.data.frame(ratio[, colnames(ratio) %in% c(samples_LPS, samples_unstim)])

ratio_m <- melt(as.matrix(ratio_LPS_unstim))
## Warning in melt(as.matrix(ratio_LPS_unstim)): The melt generic in data.table
## has been passed a matrix and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(as.matrix(ratio_LPS_unstim)). In the next version, this warning
## will become an error.
gene2check_charac = merge(ratio_m, metadata4deg, by.x = "Var2", by.y = "Sample")

# show direction of effect for male versus female
ggplot(gene2check_charac, aes(x = Stimulation, y = value, fill = Stimulation)) +
  geom_boxplot(notch = F, na.rm = T) +
  geom_jitter() +
  scale_y_continuous(limits = c(0,1)) +
  theme_bw() + facet_wrap(~Var1, scales = "free_y") +
  ggpubr::stat_compare_means(label = "p.format", label.x.npc = "centre", method = "wilcox.test")
## Warning: Removed 89 rows containing missing values (geom_point).

p <- ggplot(gene2check_charac, aes(x = Stimulation, y = value, fill = Stimulation)) + 
    geom_violin(width = 0.5) +
    geom_jitter(width = 0.25, aes(colour = Stimulation), size = 1) + 
    geom_boxplot(fill = NA, outlier.colour = NA, width = 0.25) +
    facet_wrap(~Var1) + theme_classic() 
p

Boxplot for WARS in IFNy

IFNy <- all_res$IFNy
markers = "~/Documents/MiGASti/docs/2nd_pass/WARS.xlsx"
markers = read_excel(markers, col_names = TRUE) 
markers = as.data.frame(markers)

isoform_id <- rownames(tx_ratio)
copy <- cbind (isoform_id, tx_ratio)
ratio <- merge(copy, markers, by = "isoform_id")
sample <- colnames(ratio)
rownames(ratio) = ratio$isoform_id
metadata <- read.table("~/Documents/MiGASti/Databases/metadata.txt")
#set rownames to Sample
row.names(metadata) <- metadata$Sample 
setwd("~/Documents/MiGASti/Databases")
exclude <- read.table("samples2remove.txt")
exclude <- exclude$x
ratio = ratio[, !colnames(ratio) %in% exclude] 
#Excludes the samples from filters. 
#dim(genes_tpm_filt)
metadata_filt = metadata[ !(rownames(metadata) %in% exclude), ]
length(metadata_filt)
## [1] 38
samples_IFNy = metadata_filt$Sample[metadata_filt$Stimulation %in% "IFNy"]  
samples_unstim = metadata_filt$Sample[metadata_filt$Stimulation %in% "unstim"] 

metadata4deg = metadata_filt[rownames(metadata_filt) %in% c(samples_IFNy,samples_unstim) ,]

ratio_IFNy_unstim = as.data.frame(ratio[, colnames(ratio) %in% c(samples_IFNy, samples_unstim)])

ratio_m <- melt(as.matrix(ratio_IFNy_unstim))
## Warning in melt(as.matrix(ratio_IFNy_unstim)): The melt generic in data.table
## has been passed a matrix and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(as.matrix(ratio_IFNy_unstim)). In the next version, this warning
## will become an error.
gene2check_charac = merge(ratio_m, metadata4deg, by.x = "Var2", by.y = "Sample")

# show direction of effect for male versus female
ggplot(gene2check_charac, aes(x = Stimulation, y = value, fill = Stimulation)) +
  geom_boxplot(notch = F, na.rm = T) +
  geom_jitter() +
  scale_y_continuous(limits = c(0,1)) +
  theme_bw() + facet_wrap(~Var1, scales = "free_y") +
  ggpubr::stat_compare_means(label = "p.format", label.x.npc = "centre", method = "wilcox.test")
## Warning: Removed 9 rows containing non-finite values (stat_compare_means).
## Warning: Removed 243 rows containing missing values (geom_point).

p <- ggplot(gene2check_charac, aes(x = Stimulation, y = value, fill = Stimulation)) + 
    geom_violin(width = 0.5) +
    geom_jitter(width = 0.25, aes(colour = Stimulation), size = 1) + 
    geom_boxplot(fill = NA, outlier.colour = NA, width = 0.25) +
    facet_wrap(~Var1) + theme_classic() 
p 
## Warning: Removed 9 rows containing non-finite values (stat_ydensity).
## Warning: Removed 9 rows containing non-finite values (stat_boxplot).
## Warning: Removed 9 rows containing missing values (geom_point).

Boxplot CIITA IFNy

IFNy <- all_res$IFNy
markers = "~/Documents/MiGASti/docs/2nd_pass/CIITA.xlsx"
markers = read_excel(markers, col_names = TRUE) 
markers = as.data.frame(markers)

isoform_id <- rownames(tx_ratio)
copy <- cbind (isoform_id, tx_ratio)
ratio <- merge(copy, markers, by = "isoform_id")
sample <- colnames(ratio)
rownames(ratio) = ratio$isoform_id
metadata <- read.table("~/Documents/MiGASti/Databases/metadata.txt")
#set rownames to Sample
row.names(metadata) <- metadata$Sample 
setwd("~/Documents/MiGASti/Databases")
exclude <- read.table("samples2remove.txt")
exclude <- exclude$x
ratio = ratio[, !colnames(ratio) %in% exclude] 
#Excludes the samples from filters. 
#dim(genes_tpm_filt)
metadata_filt = metadata[ !(rownames(metadata) %in% exclude), ]
length(metadata_filt)
## [1] 38
samples_IFNy = metadata_filt$Sample[metadata_filt$Stimulation %in% "IFNy"]  
samples_unstim = metadata_filt$Sample[metadata_filt$Stimulation %in% "unstim"] 

metadata4deg = metadata_filt[rownames(metadata_filt) %in% c(samples_IFNy,samples_unstim) ,]

ratio_IFNy_unstim = as.data.frame(ratio[, colnames(ratio) %in% c(samples_IFNy, samples_unstim)])

ratio_m <- melt(as.matrix(ratio_IFNy_unstim))
## Warning in melt(as.matrix(ratio_IFNy_unstim)): The melt generic in data.table
## has been passed a matrix and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(as.matrix(ratio_IFNy_unstim)). In the next version, this warning
## will become an error.
gene2check_charac = merge(ratio_m, metadata4deg, by.x = "Var2", by.y = "Sample")

# show direction of effect for male versus female
ggplot(gene2check_charac, aes(x = Stimulation, y = value, fill = Stimulation)) +
  geom_boxplot(notch = F, na.rm = T) +
  geom_jitter() +
  scale_y_continuous(limits = c(0,1)) +
  theme_bw() + facet_wrap(~Var1, scales = "free_y") +
  ggpubr::stat_compare_means(label = "p.format", label.x.npc = "centre", method = "wilcox.test")
## Warning: Removed 4 rows containing non-finite values (stat_compare_means).
## Warning: Removed 19 rows containing missing values (geom_point).

p <- ggplot(gene2check_charac, aes(x = Stimulation, y = value, fill = Stimulation)) + 
    geom_violin(width = 0.5) +
    geom_jitter(width = 0.25, aes(colour = Stimulation), size = 1) + 
    geom_boxplot(fill = NA, outlier.colour = NA, width = 0.25) +
    facet_wrap(~Var1) + theme_classic() 
p 
## Warning: Removed 4 rows containing non-finite values (stat_ydensity).
## Warning: Removed 4 rows containing non-finite values (stat_boxplot).
## Warning: Removed 4 rows containing missing values (geom_point).

LPS

empirical FDR < 0.05

length(which(all_res$LPS$empirical_FDR < 0.05))
## [1] 79

empirical FDR < 0.10

length(which(all_res$LPS$empirical_FDR < 0.10))
## [1] 119

empirical FDR < 0.15

length(which(all_res$LPS$empirical_FDR < 0.15))
## [1] 168

Table for download empirical FDR < 0.05

LPS_t <- subset(all_res$LPS, empirical_FDR < 0.05)
res_name_LPS_top = LPS_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_LPS_top)

IFNy

empirical FDR < 0.05

length(which(all_res$IFNy$empirical_FDR < 0.05))
## [1] 65

empirical FDR < 0.10

length(which(all_res$IFNy$empirical_FDR < 0.10))
## [1] 96

empirical FDR < 0.15

length(which(all_res$IFNy$empirical_FDR < 0.15))
## [1] 111

Table for download empirical FDR < 0.05

IFNy_t <- subset(all_res$IFNy, empirical_FDR < 0.05)
res_name_IFNy_top = IFNy_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_IFNy_top)

R848

empirical FDR < 0.05

length(which(all_res$R848$empirical_FDR < 0.05))
## [1] 59

empirical FDR < 0.10

length(which(all_res$R848$empirical_FDR < 0.10))
## [1] 80

empirical FDR < 0.15

length(which(all_res$R848$empirical_FDR < 0.15))
## [1] 115

Table for download empirical FDR < 0.05

R848_t <- subset(all_res$R848, empirical_FDR < 0.05)
res_name_R848_top = R848_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_R848_top)

TNFa

empirical FDR < 0.05

length(which(all_res$TNFa$empirical_FDR < 0.05))
## [1] 17

empirical FDR < 0.10

length(which(all_res$TNFa$empirical_FDR < 0.10))
## [1] 29

empirical FDR < 0.15

length(which(all_res$TNFa$empirical_FDR < 0.15))
## [1] 83

Table for download empirical FDR < 0.05

TNFa_t <- subset(all_res$TNFa, empirical_FDR < 0.05)
res_name_TNFa_top = TNFa_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_TNFa_top)

DEX

empirical FDR < 0.05

length(which(all_res$DEX$empirical_FDR < 0.05))
## [1] 2

empirical FDR < 0.10

length(which(all_res$DEX$empirical_FDR < 0.10))
## [1] 3

empirical FDR < 0.15

length(which(all_res$DEX$empirical_FDR < 0.15))
## [1] 5

Table for download empirical FDR < 0.05

DEX_t <- subset(all_res$DEX, empirical_FDR < 0.05)
res_name_DEX_top = DEX_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_DEX_top)

IL4

empirical FDR < 0.05

length(which(all_res$IL4$empirical_FDR < 0.05))
## [1] 48

empirical FDR < 0.10

length(which(all_res$IL4$empirical_FDR < 0.10))
## [1] 81

empirical FDR < 0.15

length(which(all_res$IL4$empirical_FDR < 0.15))
## [1] 147

Table for download empirical FDR < 0.05

IL4_t <- subset(all_res$IL4, empirical_FDR < 0.05)
res_name_IL4_top = IL4_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_IL4_top)

ATP

empirical FDR < 0.05

length(which(all_res$ATP$empirical_FDR < 0.05))
## [1] 509

empirical FDR < 0.10

length(which(all_res$ATP$empirical_FDR < 0.10))
## [1] 734

empirical FDR < 0.15

length(which(all_res$ATP$empirical_FDR < 0.15))
## [1] 975

Table for download empirical FDR < 0.05

ATP_t <- subset(all_res$ATP, empirical_FDR < 0.05)
res_name_ATP_top = ATP_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_ATP_top)

Age

empirical FDR < 0.05

length(which(all_res$age$empirical_FDR < 0.05))
## [1] 34

empirical FDR < 0.10

length(which(all_res$age$empirical_FDR < 0.10))
## [1] 108

empirical FDR < 0.15

length(which(all_res$age$empirical_FDR < 0.15))
## [1] 199

Table for download empirical FDR < 0.05

age_t <- subset(all_res$age, empirical_FDR < 0.05)
res_name_age_top = age_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_age_top)